/******************************************************************************/
/*                                                                            */
/*                     Basic mathematical functions                           */
/*                       for use in C programming                             */
/*                                                                            */
/*                         by C. Bach 07-2011                                 */
/*                        bach@gfz-potsdam.de                                 */
/*                                                                            */
/*                                                                            */
/******************************************************************************/

#include <stdio.h>
#include <math.h>
#include "nrutil.h"

/******************************************************************/
/*                                                                */
/*                          DMIN                                  */
/*                                                                */
/* Return minimum value of double vector                          */
/* if more than one minimum, first minimum found will be returned */
/* Input: double vector, length of vector                         */
/*                                                                */
/*    Index of first element will be 0                            */
/*                                                                */
/******************************************************************/
double dmin(long length, double *vector, long *index) {

	double minimum;
	long i;
	
	minimum = vector[0];
	*index   = 0;

	for (i=1; i<=length-1; i++) {
		if (vector[i]<minimum) {
			minimum = vector[i];
			*index   = i;
		}
	}
	
	return(minimum);
}

/******************************************/
/*                                        */
/*                 IMIN                   */
/*                                        */
/* Return minimum value of int vector     */
/* Input: int vector, length of vector    */
/******************************************/
int imin(int *vector, long length) {

	int minimum;
	int i;
	
	minimum = vector[0];

	for (i=1; i<=length-1; i++) {
		if (vector[i]<=minimum) {
			minimum = vector[i];
		}
	}
	
	return(minimum);
}

/******************************************/
/*                                        */
/*                 DMAX                   */
/*                                        */
/* Return maximum value of double vector  */
/* Input: double vector, length of vector */
/******************************************/
double dmax(double *vector, long length) {

	double maximum;
	int i;
	
	maximum = vector[0];

	for (i=1; i<=length-1; i++) {
		if (vector[i]>=maximum) {
			maximum = vector[i];
		}
	}
	
	return(maximum);
}


/******************************************/
/*                                        */
/*                 IMAX                   */
/*                                        */
/* Return maximum value of int vector     */
/* Input: int vector, length of vector    */
/******************************************/
int imax(int *vector, long length) {

	int maximum;
	int i;
	
	maximum = vector[0];

	for (i=1; i<=length-1; i++) {
		if (vector[i]>=maximum) {
			maximum = vector[i];
		}
	}
	
	return(maximum);
}

/******************************************/
/*                                        */
/*                DMEAN                   */
/*                                        */
/* Return mean value of double vector     */
/* Input: double vector, length of vector */
/******************************************/
double dmean(double *vector, long length) {

	double mean=0.0;
	int i;
	
	for (i=0; i<=length-1; i++) {
		mean += vector[i];
	}

	mean /= (double)length;
	
	return(mean);
}

/******************************************/
/*                                        */
/*                IMEAN                   */
/*                                        */
/* Return mean value of int vector        */
/* Input: int vector, length of vector    */
/******************************************/
double imean(int *vector, long length) {

	double mean;
	int i, sum=0;
	
	for (i=0; i<=length-1; i++) {
		sum += vector[i];
	}

	mean = (double)sum/(double)length;
	
	return(mean);
}

/******************************************/
/*                                        */
/*                 DSUM                   */
/*                                        */
/* Return sum of double vector            */
/* Input: double vector, length of vector */
/******************************************/
double dsum(double *vector, long length) {

	double sum=0.0;
	int i;
	
	for (i=0; i<=length-1; i++) {
		sum += vector[i];
	}
	
	return(sum);
}

/******************************************/
/*                                        */
/*                 ISUM                   */
/*                                        */
/* Return sum of integer vector           */
/* Input: int vector, length of vector    */
/******************************************/
int isum(int *vector, long length) {

	int sum=0;
	int i;
	
	for (i=0; i<=length-1; i++) {
		sum += vector[i];
	}
	
	return(sum);
}


/******************************************************/
/*                                                    */
/*                LSFIT                               */
/*                                                    */
/* linear vertical least squares fit for (x,y) data   */
/* Input: double x and y vector, number of elements   */
/* Output: constants a,b for: f(x) = a + bx           */
/******************************************************/

int lsfit(double *x, double *y, long n, double *a, double *b){

	double sumx=0, sumy=0, sumxy=0, sumx2=0;
	int i;

	for (i=0; i<=n-1; i++) {
		sumx += x[i];
		sumy += y[i];
		sumxy += x[i]*y[i];
		sumx2 += x[i]*x[i];
	}

	*a = (sumy*sumx2 - sumx*sumxy) / (n*sumx2 - sumx*sumx);
	*b = (n*sumxy - sumx*sumy) / (n*sumx2 - sumx*sumx);

	return(0);
}



/* Given an array of data[0..n-1], this routine returns:
   mean               ave
   average deviation  adev
   standard deviation sdev
   variance           var
   skewness           skew
   kurtosis           curt

   returns 1 if successful, 0 if error happened

*/

int moment(double data[], int n, double *mean, double *adev, double *sdev,
            double *var, double *skew, double *curt) {

	int j;
	double ep=0.0,s,p;

	if (n <= 1) {
		printf("\nERROR: n must be at least 2 in moment\n");
		return(0);
	}

	s=0.0;	//First pass to get the mean.

	for (j=0;j<=n-1;j++) s += data[j];

	*mean=s/n;
	*adev=(*var)=(*skew)=(*curt)=0.0;		// Second pass to get the first (absolute), sec-
	for (j=0;j<=n-1;j++) {					// ond, third, and fourth moments of the
		*adev += fabs(s=data[j]-(*mean));	// deviation from the mean.
		ep += s;
		*var += (p=s*s);
		*skew += (p *= s);
		*curt += (p *= s);
	}

	*adev /= n;
	*var=(*var-ep*ep/n)/(n-1);			// Corrected two-pass formula.
	*sdev=sqrt(*var);					// Put the pieces together according to the con-
	if (*var) {							// ventional definitions.
		*skew /= (n*(*var)*(*sdev));
		*curt=(*curt)/(n*(*var)*(*var))-3.0;
	} 
	else {
		printf("\nERROR: No skew/kurtosis when variance = 0 (in moment)\n");
		return(0);
	}
	
	return(1);
}



/******************************************/
/*                                        */
/*               DSHIFT                   */
/*                                        */
/* Shift position of two values           */
/* Input: double a, double b              */
/******************************************/
int dshift(double *a, double *b) {

	double dum;
	
	dum = *a;
	*a   = *b;
	*b   = dum;
	
	return(0);
}




/***********************************************/
/*                                             */
/*                MMULT                        */
/*                                             */
/* Matrix multiplication for symmetric matrizes*/
/* Input: double **A, double **B, double **C   */
/*        int m, int n -> no. of lines/columns */
/*                                             */
/***********************************************/
void mmult(double **A, double **B, double **C, int m, int n){
	double dum = 0.0;
	int i,j,k;

	for (i=0;i<m;i++) {
		for (j=0;j<n;j++) {
			for (k=0;k<m;k++) {
				dum += A[i][k]*B[k][j];
			}
			C[i][j] = dum;
			dum = 0.0;
		}
	}

}

/***********************************************/
/*                                             */
/*                MTRANS                       */
/*                                             */
/* Matrix transponation                        */
/* Input: double **A, double **B               */
/*        int m, int n -> no. of lines/columns */
/*                                             */
/***********************************************/
void mtrans(double **A, double **B, int m, int n) {
	int i,j;
	
	for (i=0;i<m;i++) {
		for (j=0;j<n;j++) {
			B[j][i] = A[i][j];
		}
	}

}

/***********************************************/
/*                                             */
/*             PRINTMATRIX                     */
/*                                             */
/***********************************************/
void printmatrix(double **A,int m,int n) {
	int i,j;
	
	for (i=0;i<m;i++) {
		printf("|");
		for (j=0;j<n;j++) {
			printf("%15.2lf ", A[i][j]);
		}
		printf("|\n");
	}

}


/**************************************************************/
/*                                                            */
/*                MROT                                        */
/*                                                            */
/*  Tensor Rotation around z-axis                             */
/*  degree > 0: rotate the coordinate system clockwise        */
/*  degree < 0: rotate the coordinate system counterclockwise */
/*                                                            */
/**************************************************************/
void mrot3(double **Stress, double **Result, double degree) {
	double **Rot, **Rottrans, **Mdum, deg2rad, radian;
	
	deg2rad = 0.017453292519943295;
	radian = degree*deg2rad;
	
	Rot = dmatrix(0,2,0,2);
	Rottrans = dmatrix(0,2,0,2);
	Mdum = dmatrix(0,2,0,2);

	Rot[0][0] = cos(radian);
	Rot[1][0] = -sin(radian);
	Rot[0][1] = sin(radian);
	Rot[1][1] = cos(radian);
	Rot[2][2] = 1;
	Rot[2][0] = 0.0;
	Rot[0][2] = 0.0;
	Rot[2][1] = 0.0;
	Rot[1][2] = 0.0;

	mtrans(Rot,Rottrans,3,3);
	
	mmult(Rot,Stress,Mdum,3,3);
	
	mmult(Mdum,Rottrans,Result,3,3);

	free_dmatrix(Rot,0,2,0,2);
	free_dmatrix(Rottrans,0,2,0,2);
	free_dmatrix(Mdum,0,2,0,2);
	
}


void madd(double **A, double **B, double **C, int m, int n) {
	int i,j;

	for (i=0;i<m;i++) {
		for (j=0;j<n;j++) {
			C[i][j] = A[i][j] + B[i][j];
		}
	}

}




